% matlab example of VAR
% > Econometric Toolbox > Model Selection > Specification Testing

% 1. load data
clc; clear all
load Data_USEconModel
GDP = Dataset.GDP;
M1 = Dataset.M1SL;
TB3 = Dataset.TB3MS;
Y = [GDP,M1,TB3];

%% Transforming Data for Stationarity
% plot the data for trend
figure
subplot(3,1,1)
plot(dates,Y(:,1),'r');
title('GDP')
datetick('x')
grid on
subplot(3,1,2);
plot(dates,Y(:,2),'b');
title('M1')
datetick('x')
grid on
subplot(3,1,3);
plot(dates, Y(:,3), 'k')
title('3-mo T-bill')
datetick('x')
grid on
hold off


% transfromed and plot again
Y = [diff(log(Y(:,1:2))), Y(2:end,3)]; % Transformed data
X = dates(2:end);

figure
subplot(3,1,1)
plot(X,Y(:,1),'r');
title('GDP')
datetick('x')
grid on
subplot(3,1,2);
plot(X,Y(:,2),'b');
title('M1')
datetick('x')
grid on
subplot(3,1,3);
plot(X, Y(:,3),'k'),
title('3-mo T-bill')
datetick('x')
grid on

% plot in one figure only

Y(:,1:2) = 100*Y(:,1:2);
figure
plot(X,Y(:,1),'r');
hold on
plot(X,Y(:,2),'b');
datetick('x')
grid on
plot(X,Y(:,3),'k');
legend('GDP','M1','3-mo T-bill');
hold off

%% selecting and fitting data

% === SELECTING MODEL
% make the series same length
dGDP = 100*diff(log(GDP(49:end)));
dM1 = 100*diff(log(M1(49:end)));
dT3 = diff(TB3(49:end));
Y = [dGDP dM1 dT3];

%create 4 model
dt = logical(eye(3));
VAR2diag = vgxset('ARsolve',repmat({dt},2,1),...
    'asolve',true(3,1),'Series',{'GDP','M1','3-mo T-bill'});
VAR2full = vgxset(VAR2diag,'ARsolve',[]);
VAR4diag = vgxset(VAR2diag,'nAR',4,'ARsolve',repmat({dt},4,1));
VAR4full = vgxset(VAR2full,'nAR',4);

% === Choosing Presample, Estimation, and Forecast Periods. 

Ypre = Y(1:4,:);
T = ceil(.9*size(Y,1));
Yest = Y(5:T,:);
YF = Y((T+1):end,:);
TF = size(YF,1);

% fitting
[EstSpec1,EstStdErrors1,LLF1,W1] = ...
    vgxvarx(VAR2diag,Yest,[],Ypre,'CovarType','Diagonal');
[EstSpec2,EstStdErrors2,LLF2,W2] = ...
    vgxvarx(VAR2full,Yest,[],Ypre);
[EstSpec3,EstStdErrors3,LLF3,W3] = ...
    vgxvarx(VAR4diag,Yest,[],Ypre,'CovarType','Diagonal');
[EstSpec4,EstStdErrors4,LLF4,W4] = ...
    vgxvarx(VAR4full,Yest,[],Ypre);

%% Checking Model Adequacy

[isStable1,isInvertible1] = vgxqual(EstSpec1);
[isStable2,isInvertible2] = vgxqual(EstSpec2);
[isStable3,isInvertible3] = vgxqual(EstSpec3);
[isStable4,isInvertible4] = vgxqual(EstSpec4);
[isStable1,isStable2,isStable3,isStable4]

% Likelihood Ratio Tests
% count the active number parameter
[n1,n1p] = vgxcount(EstSpec1);
[n2,n2p] = vgxcount(EstSpec2);
[n3,n3p] = vgxcount(EstSpec3);
[n4,n4p] = vgxcount(EstSpec4);
% test
reject1 = lratiotest(LLF2,LLF1,n2p - n1p)
reject3 = lratiotest(LLF4,LLF3,n4p - n3p)
reject4 = lratiotest(LLF4,LLF2,n4p - n2p)

% AIC 
AIC = aicbic([LLF1 LLF2 LLF3 LLF4],[n1p n2p n3p n4p])

% Compare forecast with forecast period data
[FY1,FYCov1] = vgxpred(EstSpec1,TF,[],Yest);
[FY2,FYCov2] = vgxpred(EstSpec2,TF,[],Yest);
[FY3,FYCov3] = vgxpred(EstSpec3,TF,[],Yest);
[FY4,FYCov4] = vgxpred(EstSpec4,TF,[],Yest);
figure
vgxplot(EstSpec2,Yest,FY2,FYCov2)
% Sum Square Error
error1 = YF - FY1;
error2 = YF - FY2;
error3 = YF - FY3;
error4 = YF - FY4;

SSerror1 = error1(:)' * error1(:);
SSerror2 = error2(:)' * error2(:);
SSerror3 = error3(:)' * error3(:);
SSerror4 = error4(:)' * error4(:);
figure
bar([SSerror1 SSerror2 SSerror3 SSerror4],.5)
ylabel('Sum of squared errors')
set(gca,'XTickLabel',...
    {'AR2 diag' 'AR2 full' 'AR4 diag' 'AR4 full'})
title('Sum of Squared Forecast Errors')

vgxdisp(EstSpec2)

%% Forecasting


%% Calculating Impulse Responses
